// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                         avtDataBinningConstructor.C                       //
// ************************************************************************* //

#include <avtDataBinningConstructor.h>

#include <vtkCell.h>
#include <vtkCellData.h>
#include <vtkDataArray.h>
#include <vtkDataSet.h>
#include <vtkPointData.h>
#include <vtkUnsignedCharArray.h>

#include <ConstructDataBinningAttributes.h>

#include <avtCallback.h>
#include <avtDataBinning.h>
#include <avtDataBinningFunctionInfo.h>
#include <avtParallel.h>
#include <avtR2Faverage.h>
#include <avtR2Fminimum.h>
#include <avtR2Fmaximum.h>
#include <avtR2Fcount.h>
#include <avtR2Fpdf.h>
#include <avtR2Fstddev.h>
#include <avtR2Fsum.h>
#include <avtR2Fvariance.h>
#include <avtR2Frms.h>
#include <avtUniformBinningScheme.h>

#include <BadIndexException.h>
#include <DebugStream.h>
#include <ExpressionException.h>
#include <PipelineException.h>
#include <ImproperUseException.h>

using     std::string;
using     std::vector;


// ****************************************************************************
//  Method: avtDataBinningConstructor constructor
//
//  Purpose:
//      Defines the destructor.  Note: this should not be inlined in the header
//      because it causes problems for certain compilers.
//
//  Programmer: Hank Childs
//  Creation:   February 12, 2006
//
// ****************************************************************************

avtDataBinningConstructor::avtDataBinningConstructor()
{
    ;
}


// ****************************************************************************
//  Method: avtDataBinningConstructor destructor
//
//  Purpose:
//      Defines the destructor.  Note: this should not be inlined in the header
//      because it causes problems for certain compilers.
//
//  Programmer: Hank Childs
//  Creation:   February 12, 2006
//
// ****************************************************************************

avtDataBinningConstructor::~avtDataBinningConstructor()
{
    ;
}


// ****************************************************************************
//  Method: avtDataBinningConstructor::ConstructDataBinning
//
//  Purpose:
//      Constructs the actual data binning.
//
//  Programmer: Hank Childs
//  Creation:   February 12, 2006
//
//  Modifications:
//
//    Hank Childs, Sat Feb 25 15:01:56 PST 2006
//    Add support for time.
//
//    Hank Childs, Mon May 15 14:25:28 PDT 2006
//    Fix memory leak.
//
//    Hank Childs, Fri Jun  9 13:44:39 PDT 2006
//    Fix problem with mixed centering. Also add "default" to switch statement.
//
//    Hank Childs, Mon Aug 20 10:07:01 PDT 2007
//    Detect and report errors that happen upstream when getting the data ready
//    to calculate the DDF.
//
//    Hank Childs, Mon Dec 10 16:54:53 PST 2007
//    Make sure, for unstructured grids, we are only using nodes that are
//    referenced by cells.
//
//    Sean Ahern, Thu Jan 10 16:01:06 EST 2008
//    Added a "count" statistical operator.
//
//    Cyrus Harrison, Tue Aug 10 10:48:17 PDT 2010
//    Added a "rms" operator.
//
//    Hank Childs, Sat Aug 21 14:05:14 PDT 2010
//    Better support for mixed centering, r2f's that don't require a variable.
//
//    Hank Childs, Tue Sep 13 14:17:34 PDT 2011
//    Add support for data binning's based on spatial coordinates.
//
//    Hank Childs, Thu Sep 22 08:39:24 PDT 2011
//    Fix crash when processor has no data.
//
//    Hank Childs, Mon Jul 16 16:19:51 PDT 2012
//    Fix handling of mixed centering when the statistical reduction does not 
//    involve a variable.
//
//    Cyrus Harrison, Tue Feb 23 10:41:21 PST 2021
//    Add proper ghost zone handling (skip ghosts during binning)
//
// ****************************************************************************

class ValueRetriever
{
  public:
    virtual double GetValue(int) = 0;
    virtual ~ValueRetriever() { ; };
};

class VariableValueRetriever : public ValueRetriever
{
  public:
                   VariableValueRetriever(vtkDataArray *a) { arr = a; };
    double         GetValue(int v) { return arr->GetTuple1(v); };
  private:
    vtkDataArray  *arr;
};

class CoordinateValueRetriever : public ValueRetriever
{
  public:
                   CoordinateValueRetriever(vtkDataSet *m, int c) { mesh = m; coord = c; };
    double         GetValue(int v) { double p[3]; mesh->GetPoint(v, p);  return p[coord]; };

  private:
    vtkDataSet    *mesh;
    int            coord;
};


avtDataBinning *
avtDataBinningConstructor::ConstructDataBinning(
                                         ConstructDataBinningAttributes *atts,
                                         avtContract_p spec, bool mustReExecute)
{
    //
    // Certify the attributes.
    //
    if (atts->GetBinningScheme() != ConstructDataBinningAttributes::Uniform)
    {
        debug1 << "Only uniform binning is supported" << endl;
        EXCEPTION1(ImproperUseException, "Only uniform binning is supported");
        return NULL;
    }
    if (atts->ReductionRequiresVariable() && atts->GetVarForReductionOperator() == "")
    {
        debug1 << "The reduction operator requires a variable." << endl;
        EXCEPTION1(ImproperUseException, "The reduction operator requires a variable.");
        return NULL;
    }
    for (size_t i = 0 ; i < atts->GetVarnames().size() ; i++)
    {
        if (atts->GetVarnames()[i] == "")
        {
            debug1 << "Some variable names are not valid." << endl;
            EXCEPTION1(ImproperUseException, 
                       "Some variable names are not valid.");
            return NULL;
        }
    }
    if ((atts->GetBinType().size() != atts->GetNumBins().size()) &&
        (atts->GetBinType().size() != 0))
    {
        debug1 << "Bin types were specified, but the wrong number were "
               << "specified." << endl;
        EXCEPTION1(ImproperUseException,
                   "If you specify bin types, the number of bin types must "
                   "match the number of variables.");
        return NULL;
    }

    // For legacy reasons, binType might be a zero-length array.  In that 
    // case, the binType should be assumed to be all "Variable".  Make an
    // array following this assumption, so coding below can be a bit easier.
    int *binType = new int[atts->GetNumBins().size()];
    for (size_t i = 0 ; i < atts->GetNumBins().size() ; i++)
    {
        if (atts->GetBinType().size() == 0)
            binType[i] = ConstructDataBinningAttributes::Variable;
        else
            binType[i] = atts->GetBinType()[i];
    }

    if (atts->GetNumBins().size() != atts->GetVarnames().size())
    {
        debug1 << "There must be the same number of samples "
                                  "as there are variable names." << endl;
        EXCEPTION1(ImproperUseException,
                   "There must be the same number of samples "
                   "as there are variable names.");
        return NULL;
    }
    if (atts->GetBinBoundaries().size() != 2*atts->GetVarnames().size())
    {
        debug1 << "There must be a minimum and a maximum for "
                  "every variable.  The range is the wrong size." << endl;
        EXCEPTION1(ImproperUseException,
                   "There must be a minimum and a maximum for "
                   "every variable.  The range is the wrong size.");
        return NULL;
    }
    for (size_t i = 0 ; i < atts->GetVarnames().size() ; i++)
    {
        if (atts->GetBinBoundaries()[2*i] >= atts->GetBinBoundaries()[2*i+1])
        {
            debug1 << "The minimum for the range of a variable"
                                      "must be less than the maximum" << endl;
            EXCEPTION1(ImproperUseException,
                       "The minimum for the range of a variable"
                       "must be less than the maximum");
            return NULL;
        }
    }
    for (size_t i = 0 ; i < atts->GetVarnames().size() ; i++)
    {
        if (atts->GetNumBins()[i] <= 0)
        {
            debug1 << "Each variable must have at least one bin" << endl;
            EXCEPTION1(ImproperUseException,
                       "Each variable must have at least one bin");
            return NULL;
        }
    }
   

    //
    // Set up the binning scheme and the function info based on the input
    // arguments.
    //
    int nvars = atts->GetVarnames().size();
    float *minmax = new float[2*nvars];
    for (size_t i = 0 ; i < (size_t)2*nvars ; i++)
        minmax[i] = atts->GetBinBoundaries()[i];
    int *nvals = new int[nvars];
    for (size_t i = 0 ; i < (size_t)nvars ; i++)
        nvals[i] = atts->GetNumBins()[i];
    avtBinningScheme *bs = new avtUniformBinningScheme(nvars, minmax, nvals);
    std::vector<std::string> varnames = atts->GetVarnames();
    avtDataBinningFunctionInfo *info = new avtDataBinningFunctionInfo(bs, 
             varnames, atts->GetBinType(), atts->GetVarForReductionOperator());

    //
    // Now create the "relation-to-function operator".
    //
    int nBins = bs->GetNumberOfBins();
    bs->SetOutOfBoundsBehavior(atts->GetOutOfBoundsBehavior());
    avtR2Foperator *R2Foperator = NULL;
    switch (atts->GetReductionOperator())
    {
      case ConstructDataBinningAttributes::Average:
        R2Foperator = new avtR2Faverage(nBins, atts->GetUndefinedValue());
        break;
      case ConstructDataBinningAttributes::Minimum:
        R2Foperator = new avtR2Fminimum(nBins, atts->GetUndefinedValue());
        break;
      case ConstructDataBinningAttributes::Maximum:
        R2Foperator = new avtR2Fmaximum(nBins, atts->GetUndefinedValue());
        break;
      case ConstructDataBinningAttributes::Count:
        R2Foperator = new avtR2Fcount(nBins);
        break;
      case ConstructDataBinningAttributes::PDF:
        R2Foperator = new avtR2Fpdf(nBins);
        break;
      case ConstructDataBinningAttributes::StandardDeviation:
        R2Foperator = new avtR2Fstddev(nBins, atts->GetUndefinedValue());
        break;
      case ConstructDataBinningAttributes::Sum:
        R2Foperator = new avtR2Fsum(nBins);
        break;
      case ConstructDataBinningAttributes::Variance:
        R2Foperator = new avtR2Fvariance(nBins, atts->GetUndefinedValue());
        break;
      case ConstructDataBinningAttributes::RMS:
        R2Foperator = new avtR2Frms(nBins, atts->GetUndefinedValue());
        break;
      default:
        // Caught in logic below
        break;
    }
    if (R2Foperator == NULL)
    {
        EXCEPTION2(ExpressionException, atts->GetName(),
                   "You have requested an operator that has not been implemented yet.");
    }

    //
    // Make the input pipeline execute again to get the variables we need.
    //
    int timeStart = spec->GetDataRequest()->GetTimestep();
    int timeEnd = timeStart+1;
    int timeStride = 1;
    if (atts->GetOverTime())
    {
        timeStart  = atts->GetTimeStart();
        timeEnd    = atts->GetTimeEnd();
        timeStride = atts->GetTimeStride();
    }

    //
    // Now iterate some number of passes over the data.  This will probably
    // just be one pass, but could be multiple passes in the case of operators
    // like "standard deviation".  On the final pass, the operator will return
    // the actual data.  On previous passes, it will return NULL.
    //
    int nPasses = R2Foperator->GetNumberOfPasses();
    float *vals = NULL;
    for (int pass = 0 ; pass < nPasses ; pass++)
    {
        for (int time = timeStart ; time < timeEnd ; time += timeStride)
        {
            avtContract_p spec2 = new avtContract(spec);
            avtDataRequest_p dataRequest = spec2->GetDataRequest();
            if (atts->GetVarForReductionOperator() != "")
                dataRequest->AddSecondaryVariable(atts->GetVarForReductionOperator().c_str());
            for (size_t i = 0 ; i < atts->GetVarnames().size() ; i++)
                if (binType[i] == ConstructDataBinningAttributes::Variable)
                    dataRequest->AddSecondaryVariable(atts->GetVarnames()[i].c_str());
            dataRequest->SetTimestep(time);
            if (mustReExecute)
                GetInput()->Update(spec2);
            if (GetInput()->GetInfo().GetValidity().HasErrorOccurred())
            {
                char msg[1024];
                snprintf(msg, 1024, "The data binning could not be calculated because "
                                    "an error occurred when generating "
                                    "the data to construct it from.  The error"
                                    " was \"%s\".",
                GetInput()->GetInfo().GetValidity().GetErrorMessage().c_str());
                EXCEPTION1(ImproperUseException, msg);
            }

            //
            // Get all of the input data.
            //
            avtDataTree_p tree = GetInputDataTree();
            int nLeaves;
            vtkDataSet **leaves = tree->GetAllLeaves(nLeaves);
        
            bool coDomIsNodal = true;
            bool *isNodal = new bool[nvars];
            const char   *codomain_varname = info->GetCodomainName().c_str();
            bool hasError = false;
            bool mixedCentering = false;
            if (nLeaves > 0)
            {
                for (size_t k = 0 ; k < (size_t)nvars ; k++)
                {
                    avtDataBinningFunctionInfo::BinBasedOn bbo = 
                                                    info->GetBinBasedOnType(k);
                    if (bbo == avtDataBinningFunctionInfo::VARIABLE)
                    {
                        const char *varname=info->GetDomainTupleName(k).c_str();
                        if (leaves[0]->GetPointData()->GetArray(varname)!=NULL)
                            isNodal[k] = true;
                        else if (leaves[0]->GetCellData()->GetArray(varname) 
                                                                        != NULL)
                            isNodal[k] = false;
                        else
                            hasError = true;
                    }
                    else
                    {
                        // Coordinate, so nodal
                        isNodal[k] = true;
                    }
                }
    
                if (atts->ReductionRequiresVariable())
                {
                    if (leaves[0]->GetPointData()->GetArray(codomain_varname)
                                                                       != NULL)
                        coDomIsNodal = true;
                    else if (leaves[0]->GetCellData()->GetArray(codomain_varname) 
                                                                       != NULL)
                        coDomIsNodal = false;
                    else
                        hasError = true;

                    for (size_t k = 0 ; k < (size_t)nvars ; k++)
                        if (isNodal[k] != coDomIsNodal)
                            mixedCentering = true;
                }
                else
                {
                    // still should check for mixed centering
                    for (size_t k = 1 ; k < (size_t)nvars ; k++)
                        if (isNodal[k] != isNodal[0])
                            mixedCentering = true;
                }
            }

            int val = UnifyMaximumValue((hasError ? 1 : 0));
            if (val > 0)
            {
                // Free the memory from the GetAllLeaves function call.
                delete [] leaves;
                delete [] isNodal;

                debug1 << "Could not create data binning because either a "
                                  "variable could not be located, or because"
                                  " the centering of the variables do not "
                                  "match." << endl;
                EXCEPTION1(ImproperUseException, 
                                  "Could not create data binning because either a "
                                  "variable could not be located, or because"
                                  " the centering of the variables do not "
                                  "match.");
                return NULL;
            }
        
            ValueRetriever **val_ret = new ValueRetriever*[nvars];
            for (size_t k = 0 ; k < (size_t)nvars ; k++)
                 val_ret[k] = NULL;
            float        *args = new float[nvars];
            for (size_t j = 0 ; j < (size_t)nLeaves ; j++)
            {
                vtkDataArray *codomain = (coDomIsNodal 
                        ? leaves[j]->GetPointData()->GetArray(codomain_varname)
                        : leaves[j]->GetCellData()->GetArray(codomain_varname));
                for (size_t k = 0 ; k < (size_t)nvars ; k++)
                {
                    avtDataBinningFunctionInfo::BinBasedOn bbo = 
                                                    info->GetBinBasedOnType(k);
                    if (bbo == avtDataBinningFunctionInfo::VARIABLE)
                    {
                        const char *varname = info->GetDomainTupleName(k).c_str();
                        vtkDataArray *arr = (isNodal[k]
                                    ? leaves[j]->GetPointData()->GetArray(varname)
                                    : leaves[j]->GetCellData()->GetArray(varname));
                        val_ret[k] = new VariableValueRetriever(arr);
                    }
                    else if (bbo == avtDataBinningFunctionInfo::X)
                        val_ret[k] = new CoordinateValueRetriever(leaves[j], 0);
                    else if (bbo == avtDataBinningFunctionInfo::Y)
                        val_ret[k] = new CoordinateValueRetriever(leaves[j], 1);
                    else if (bbo == avtDataBinningFunctionInfo::Z)
                        val_ret[k] = new CoordinateValueRetriever(leaves[j], 2);
                }
                int nvals;
                bool doCells = true;
                if (mixedCentering)
                    doCells = true;
                else if (atts->ReductionRequiresVariable())
                {
                    if (!coDomIsNodal)
                        doCells = true;
                    else
                        doCells = false;
                }
                else
                {
                    if (! isNodal[0])
                        doCells = true;
                    else
                        doCells = false;
                }

                nvals = (doCells ? leaves[j]->GetNumberOfCells()
                                 : leaves[j]->GetNumberOfPoints());

                vector<bool> useValue;
                if (!doCells)
                {
                    if (leaves[j]->GetDataObjectType() == VTK_UNSTRUCTURED_GRID ||
                        leaves[j]->GetDataObjectType() == VTK_POLY_DATA)
                    {
                        useValue.resize(nvals, false);
                        size_t ncells = leaves[j]->GetNumberOfCells();
                        for (size_t l = 0 ; l < ncells ; l++)
                        {
                            vtkCell *cell = leaves[j]->GetCell(l);
                            for (int p = 0 ; p < cell->GetNumberOfPoints() ; p++)
                            {
                                vtkIdType ptId = cell->GetPointId(p);
                                useValue[ptId] = true;
                            }
                        }
                    }
                }

                // for cells case: 
                // see if we have a ghost zone labels array, if so
                // we want to use it to exclude ghosts when constructing
                // the binning
                vtkUnsignedCharArray *gzones = NULL;
                if(doCells)
                {
                    if(leaves[j]->GetCellData()->HasArray("avtGhostZones"))
                    {
                        debug5 << "DDF input dataset " << j 
                               << " has ghost zones" << endl;
                        gzones = (vtkUnsignedCharArray *) leaves[j]->GetCellData()->GetArray("avtGhostZones");
                    }
                    else
                    {
                        debug5 << "DDF input dataset " << j 
                               << " does not have ghost zones" << endl;
                    }
                }

                for (size_t l = 0 ; l < (size_t)nvals ; l++)
                {
                    if (useValue.size() > 0)
                        if (!useValue[l])
                            continue;

                    // skip if this zonal operation and the current zone is ghost
                    if(gzones != NULL)
                    {
                        if(gzones->GetValue(l) != 0)
                        {
                            continue;
                        }
                    }

                    if (!mixedCentering)
                    {
                        for (size_t k = 0 ; k < (size_t)nvars ; k++)
                            args[k] = val_ret[k]->GetValue(l);
                        int binId = bs->GetBinId(args);
                        float cval = 0.;
                        if (atts->ReductionRequiresVariable())
                            cval = codomain->GetTuple1(l);
                        if (binId >= 0)
                            R2Foperator->AddData(binId, cval);
                    }
                    else
                    {
                        vtkCell *cell = leaves[j]->GetCell(l);
                        for (int p = 0 ; p < cell->GetNumberOfPoints() ; p++)
                        {
                            vtkIdType ptId = cell->GetPointId(p);
                            for (size_t k = 0 ; k < (size_t)nvars ; k++)
                                if (isNodal[k])
                                    args[k] = val_ret[k]->GetValue(ptId);
                                else
                                    args[k] = val_ret[k]->GetValue(l);
                            int binId = bs->GetBinId(args);
                            float cval = 0;
                            if (atts->ReductionRequiresVariable())
                                 cval =  (coDomIsNodal
                                                 ? codomain->GetTuple1(ptId)
                                                 : codomain->GetTuple1(l));
                            if (binId >= 0)
                                R2Foperator->AddData(binId,cval);
                        }
                    }
                }
            }
            for (size_t j = 0 ; j < (size_t)nvars ; j++)
                if (val_ret[j] != NULL)
                    delete val_ret[j];
            delete [] val_ret;
            delete [] args;
            delete [] leaves;
            delete [] isNodal;
        }

        vals = R2Foperator->FinalizePass(pass);
    }

    //
    // We can now construct the actual data binning.
    //
    avtDataBinning *rv = new avtDataBinning(info, vals);

    //
    // Clean up memory.
    //
    delete R2Foperator;
    delete [] minmax;
    delete [] nvals;
    delete [] binType;
    // Don't delete binning scheme (bs) function info (info), or "vals"
    // since the data binning now owns them.

    return rv;
}


